Rare fluctuations and large-scale circulation cessations in turbulent convection 
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In turbulent Rayleigh-Benard convection, a large-scale circulation (LSC) develops in a nearly 
vertical plane, and is maintained by rising and falling plumes detaching from the unstable thermal 
boundary layers. Rare but large fluctuations in the LSC amplitude can lead to extinction of the 
LSC (a cessation event), followed by the re-emergence of another LSC with a different (random) 
azimuthal orientation. We extend previous models of the LSC dynamics to include momentum and 
thermal diffusion in the azimuthal plane, and calculate the tails of the probability distributions of 
both the amplitude and azimuthal angle. Our analytical results are in very good agreement with 
experimental data. 
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When a fluid is heated from below in the presence of a 
gravitational field, the static state with thermal conduc- 
tion can become unstable towards a succession of insta- 
bilites, ultimately leading to turbulence if the buoyancy- 
induced driving force is sufficiently greater than the vis- 
cous drag and diffusion of heat. This balance is quantified 
by the Rayleigh number Ra = a.^^gt^TlP' jvn, where a.^ is 
the isobaric thermal expansion coefficient, g is the grav- 
ity field, AT is the temperature gap between bottom and 
top layers, L is the height of the fluid container, k is the 
thermal diffusivity and v is the kinematic viscosity. For 
large Ra, thermal boundary layers become unstable by 
emitting hot (on the bottom) or cold (on the top) plumes 
which, due to buoyancy, migrate upwards (hot) or down- 
wards (cold) 1-3 . In addition to their vertical motion, 
plumes drift along the top and bottom boundaries in op- 
posite directions, contributing to a large-scale circulation 
(LSC) flowing in a nearly vertical plane, which spans the 
diameter of the container. The horizontal velocity of the 
plumes oscillates rapidly compared to the reorientation 
dynamics of the large-scale circulation [^{B], which, in a 
cylindrical geometry, undergoes both rotational diffusion 
and orientational jumps following irregular cessation of 
the entire flow [H [7] . Such laboratory experiments pro- 
vide a well-controlled setting in which to study the sta- 
tistical properties of cessation, reversal and reorientation 
events similar to those that occur in many flows of prac- 
tical significance, including atmospheric [8] and oceanic 
circulation [9], the dynamo driving planetary magnetic 
fields Hn], and in the cores of stars llj. 

In order to interpret high quality data on the statistics 
of cessations and azimuthal rotation, a nonlinear stochas- 
tic model that retains physically relevant aspects of the 
Navier-Stokes equations was developed and shown to re- 
produce many aspects of the statistics of the azimuthal 
dynamics and the temperature fluctuations in the LSC 
plane [12l [13] . The stochastic variables in the model are 
the amplitude of azimuthal temperature variations, b, in- 
duced by the LSC and the azimuthal orientation angle. 



(?o, of the nearly vertical LSC plane. Although the model 
predictions are in good agreement with the experimen- 
tal results for typical fluctuations of the system P3] , the 
model does not account quantitatively for the rare large 
fluctuations responsible for the cessation statistics and 
for the broad-tail probability distribution function (PDF) 
of the azimuthal velocities. 

The purpose of this Letter is to extend the stochastic 
model to capture the tail of the PDFs of the tempera- 
ture amplitudes and azimuthal velocities. We make three 
contributions here. Firstly, we show that the equation 
for the amplitude b needs to explicitly include a constant 
term, known to scale as Ra^^"*. Such a term was already 
proposed in Ref. [13j as arising from boundary layer ther- 
mal diffusion, but its significance for the asymptotics of 
the PDF had not been emphasized. Secondly, we show 
that the description of the azimuthal velocities needs to 
include viscous diffusion in the boundary layer near the 
wall. Such a term is generally small compared to the 
other terms in the equation of motion for but be- 
comes the dominant contribution when the amplitude b 
is small, as in a cessation event. Thirdly, we compute 
the PDFs for both b and ^q, predicting respectively an 
exponential dependence at small b and a power law of 
—4 for the large angular velocity asymptotics of ^o- A 
careful analysis of the experimental data is in very good 
agreement with these predictions. 

Evolution equation for the LSC amplitude:- For com- 
pleteness, we briefly summarize the derivation of the 
physical model for LSC fluctuations, largely following 
Refs. [m [13], but with minor differences noted below. 
The LSC amplitude evolution is derived from the equa- 
tion satisfied by the velocity component in the LSC plane, 
M0, where only the buoyancy and diffusion terms are re- 
tained. Here, <j) is the angle in the vertical circulation 
plane of the LSC. The turbulent advection term is dis- 
carded on the basis that the convection due to azimuthal 
motion is small relative to the other terms and the self- 
advection is replaced by random fiuctuations. A spatial 
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average in a direction perpendicular to the main axis of 
the cylinder (radial average) is performed. The buoy- 
ancy term acts everywhere in the LSC plane, hence the 
average keeps the same form. On the other hand, mo- 
mentum diffusion is assumed to dominate only in the 
viscous boundary layer, so that the average of this term 
gives a prefactor X/L, where A is the viscous boundary 
layer thickness. The viscous layer thickness can be esti- 
mated on dimensional grounds as the length scale where 
the convective forces balance out the diffusive forces, giv- 
ing V'^jL ~ vU/\^, where U{t) is the maximum speed 
just within the viscous boundary layer, and thus is an 
estimate for the typical turnover velocity of an eddy 
spanning the LS'C-plane. Hence A(t) ^ \fvLj\J and 
the spatial average of the diffusion term is estimated as 
{vXl'^u^) ~ -vU/{LX) ~ -z/1/2[/3/2/l3/2. Furthermore, 
we assume that the amplitude of azimuthal temperature 
variation, 5, is proportional to the large-scale typical ve- 
locity of thermal convection rolls U , in the approximation 
that momentum acceleration is due to buoyancy forces; 
this proportionality argument is different from one used 
in Ref. |13| . where buoyancy is balanced against diffu- 
sion. Finally, a delta-correlated Gaussian stochastic forc- 
ing fs{t) with amplitude Ds is included to simulate the 
effect of turbulent fluctuations. As noted in Ref. 13], 
the resulting equation incorrectly accounts for the small 
5 behavior, where the thermal boundary layer cannot be 
neglected. Thermal diffusion leads to a constant driving 
term 6 = A empirically |7j found to scale as Ra^^''. This 
has the effect of driving the system back to the vicinity 
of (5 where the mean LSC amplitude is denoted by 

5q w ATcrRe'^^^/Ra and a = v/k is the Prandtl number. 
By rescaling time t — > t/rs, where ts ~ L"^ /{lyKe^^^} is 
the typical turnover time, and defining a dimensionless 
amplitude ^ = S/Sq, we arrive at the following Langevin 
equation for the LSC fluctuations ^ 



(1) 



with A = Ats/5[). Here the rescaled (dimensionless) 
diffusion coefficient is Ds = Dsts/Sq, representing the 
amplitude of the scaled noise f({t). We have included 
numerical prefactors a,/3 = 0{1) to account for the geo- 
metric coefficients from the spatial volume averaging pro- 
cedure. These constants will be determined below by de- 
manding that the maximum and the width of the PDF 
are consistent with experimental results. 

Evolution of the azimuthal velocity:- The equation for 
the horizontal motion is obtained from the Navier-Stokes 
equation for the azimuthal velocity, ug ~ L9, by retain- 
ing the advection and momentum diffusion terms. Previ- 
ously dH [T3] , the viscous drag term was neglected on 
the basis that it is typically small. This approxima- 
tion is valid in the regime of a well-defined LSC, but 
breaks down near cessations, since the momentum trans- 
port from the LSC also becomes very small. The vis- 
cous drag is dominant in the viscous boundary layer, so 



that a spatial average along an arbitrary direction in the 
horizontal plane gives (i/V'^Bq) ~ —v9[)/{LXe). The vis- 
cous boundary layer thickness Xg is estimated from bal- 
ancing the advection force with the momentum diffusion 
force U9o/L ~ vOq/X^\ together with the proportionality 
U/ts ^ ag6, we find that Xg ~ ^/ vLSq/ {Uq5). In addi- 
tion to these deterministic forces, the self- advection term 
is mimicked by a delta-correlated Gaussian noise fg(t) 
with amplitude Dg. Rescaling time by the typical time 
Tg K, Li^ /[uRe) for crossing a boundary layer of thickness 
Xg and 5 hy 5q, the equation of motion for the azimuthal 
fluctuations is 



= - "1^ 
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where the rescaled diffusion coefficient is Dg — DgTg, 
and ai,/3i = 0(1) account for geometrical factors due to 
volume averaging. From the definition of the timescales, 
Tg/rs ~ Re~^^^ <C 1 and hence the viscous drag term 
becomes important when ^ {/3i/ai)^Ke~^ , i.e. near 
cessations. 

Probability distribution for d:- Since the Langevin 
equation for is decoupled from that of 9, we first analyze 
Eq. ([T]) separately. In order to obtain the stationary PDF 
P(^) at long times, we use the equivalent Fokker-Planck 
equation of Eq. ([!]) . It reads [H] 
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The stationary solution of this equation is 



P{0 = CeM-2ViO/Dsl 
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Note that Eq. Q predicts that logP(^ < 1) cx ^ as ob- 
served in experiment. Denoting the logarithmic deriva- 
tive of the experimental PDF at small ^ by B, using Q 
we find that A = BDsl2. Here B and Ds are the tun- 
ing parameters of the theory and will be extracted from 
experimental data. 

We now determine the constants a and /? by requiring 
that the PDF has a maximum at ^ = 1 and width equal 
to \/Ds, and fix the constant C by normalizing P(^) in 
its Gaussian regime close to ^ = 1. Expanding P(f ) Q 
in the vicinity of ^ = 1 up to second order, we find that 
A + a — j3 — for the maximum to be at ^ = 1, and 
(3/2)/3 — a ~ 1/2 for the variance to be Ds- This yields 
a = 1 - 3i and /3 = 1 - 2i. With A = BDs/2, the final 
normalized result for the PDF reads as 
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FIG. 1. (Color online) Theoretical (solid line) and experimen- 
tal (triangles) PDFs for the normalized amplitude ^ versus 

(^ — 1)/(J (cr — \/ Ds), for different Ra numbers, with (a)-(h) 
for the medium sample and (i)-(p) for the large sample. The 
Ra numbers are: 3.78 ■ 10* (a), 8.16 • 10* (b), 1.1 • lO'' (c), 
2.3-10^ (d), 4.5-10^ (e), 7.9-10^ (f), 1.02-10^° (g), 1.51-10^° 
(h), 4.75 • 10** (i), 7.16 • 10^ (j), 1.22 • 10'° (k), 2.43 • 10^" (1), 
4.71-10^° (m), 5.68-10^° (n), 7.51-10^° (o) and 1.04-10" (p). 
In each subfigure the parameters Ds and B were computed 
by fitting the left tail of the experimental PDF to Eq. ([7|. 
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FIG. 2. (Color online) Cessation frequency (per day) as a 
function of Ra for the medium (a) and large (b) samples. Ex- 
perimental results (triangles) [T3] are compared to theoretical 
prediction (10 1. The latter are bound within the two solid 
lines. The experimental cessation was defined to occur when 
S/5o < ^min, that is, we have averaged over the time intervals 
between events wh ere the system undergone cessation with 
C < ^min, see Eq. (10 1. Here, ^min = 0.15 for the medium 
sample, and ^min ~ 0.2 for the large sample. Similar results 
were obtained for thresholds of £,min between 0.15 and 0.3. 



The only free parameters in this result are B and Dg- 
These are estimated from the ^ ^ 1 asymptote of ([6]): 



P(e « 1) ~ i27rDs) 



-l/2g-3B/10-l/(5Da) ^B^ ^ 
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By fitting the logarithm of the experimental PDFs to a 
line from the logarithm of Eq. ([t]) , we extract B and Ds 
for each experimental data set. In Fig. [T] we show compar- 
isons of experimental results and PDF ^ using the pa- 
rameters B and Ds extracted from experimental results; 
very good agreement is evident for a wide range of Ra 
numbers. In this figure and henceforth the medium and 
large samples refer to cylindrical containers with heights 
24.76cm and 50.61cm and aspect ratio of ~ 1 [T5] . 

Having calculated the complete PDF of the LSC ampli- 
tudes, we now extract the cessation frequency. The lat- 
ter can be found by analyzing the following first-passage 
problem: starting from the vicinity of the fixed point 
^ = 1 what is the mean time it takes to reach the vicinity 
of ^ = ^0 ^ Ij where ^ 1 is the amplitude which de- 
fines the experimental cessation threshold? Writing down 
the backward Fokker-Planck equation Jl] , the mean time 
to cessation (MTC) is given by 

/■^ dll P°° ^b(z) 2[v(^)-v(i;o)l 

T(C,?o) = 2/-^/ mdz; ^P{z)^e , 
J^o V{v) Jy Ds 

(8) 

where ^ ~ 1 is the effective initial condition, and the 



potential satisfies ([5|) with A ~ BDs/2. Using the small- 
ness of Ds (typically ranging between 10~^-10^^), we 
can evaluate the inner integral by using the saddle point 
approximation. By doing so, we arrive at a result inde- 
pendent of y, which permits the evaluation of the outer 
integral using a Taylor expansion of the integrand about 
y = £,0 [H]. This procedure leads to the final result for 
the MTC Tc(^o) to reach a point < 1 (see also [13]): 



\V'{^o)\^Ds 



(9) 



where we have multiplied the result by ts to present the 
time in physical units, and used the fact that V"{1) = 
1/2. Given a threshold for cessation as is done ex- 
perimentally, in order to mimic the binning procedure 
of the experimental data, we have to average over in 
Eq. ([9| from to ^min, which yields 



(10) 



To obtain a theoretical prediction for ujc as a function of 
Ra, we use the extracted values oi B, Ds and ts from 
the experimental data [13]. By doing so, we can plot the 
theoretical prediction for oj^ as a function of Ra, as shown 
in Fig.[2j Here, the theoretical predictions ( 10 ) agree well 
with the experimental data |13j , for both the medium and 
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FIG. 3. (Color online) PDF P{Ae) averaged over all ex- 
perimental PDFs with Ra numbers ranging from 10^-10" 
as a function of AS/o"*; a, is a rescaling factor so that the 
medium (triangles) and large (squares) sample PDFs coincide 
at their right tail. The experimental data shown here was 
adaptively binned using the data threshold method [TB] with 
threshold 0.001 max[P( AS)]. Using the least-square method 
in the regime 4 < AS/cr* < 20, the medium and large sam- 
ples were found to scale with a power law of —4.34 ± 0.02 
and — 4.27± 0.02, respectively, closely fitting the theoretical 
prediction of Eq. (|13||. The solid line is a power law of —4.3. 



large samples. The error bars in the experimental results 
originate from the binning method, while the errors in 
the theoretical curves come from the uncertainties in the 
extracted values of B and Ds- 

Probability distribution for 9:- Now we turn to the cal- 
culation of P{0). As the Langevin equation for 6 [Eq. ^] 
depends on ^, for a given ^ we can first determine the 
steady state conditional PDF P{9q\£_). Using (121), we find 



P{0o\O = 



1 
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Given the PDF P{£_) from Eq. we then determine 
the complete PDF P{0o) of the azimuthal velocity by 
the following relation 



(12) 



which is valid when the relaxation timescale of 6 is much 
faster than that of 6, namely ts ^ tq. That is, Eq. (12) 



holds when the conditional PDF P(0olO equilibrates 
much faster than the typical timescale of change of ^. 



Integral ( 12 ) can be evaluated in the Gaussian regime 
of the PDF, where ^ — 1, which yields the statistics of 
reorientations due to rotations of the LSC plane. In this 
case, to leading order one can simply put ^ = 1 in Eq. ^ 
which gives P{0o) ^ e~°'^^<^^^i> . By comparing it with the 
experiments, ai = 1 in agreement with Ref. |13| . 

The cessation events correspond to the right hand tail 
of the PDF (12 1. Indeed, when the system undergoes 



cessation and ^ <C 1, the integrand is dominated by 



Therefore, the right hand tail of the 



PDF given by Eq. ( 12 ) satisfies 



P{Oo) - 0, 
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This power-law prediction for the tail of P{0o) is con- 
sistent with the earliest analysis of the experimental 
data [3 [12] (reporting an exponent of —3.8), but differs 
from the numerical calculations presented in [12J which 
obtained a power-law with an exponent of —2. The dif- 
ference arises as our equation of motion includes momen- 
tum diffusion that allows accurately accounting for the 
tailofP(^o)- In Fig. [s] we plot the experimental PDFs for 
A9 ^ 9, and show that they do indeed exhibit a power- 
law behavior at the tails with exponent of approximately 
—4.3 in very good agreement with our prediction. 
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